#! /bin/gawk -f

# Script to convert a Biosym mdf file into an Aqua PseudoLIB file,
# replacing residue numbers by the wildcard symbol, and removing
# multiple pseudo atom definitions.
#
# arguments: 
#    type (biosym|disgeo|diana|xplor|aqua) 
#                NOTE: code for options other than aqua has not been updated....
#                see mdf2AtomLIB
#    pdb (on|off) (default: pdb=on)
#
# mdf format:
# @list pseudoatom  PCOL:LEU_11:HBX A
# LEU_11:HB1 LEU_11:HB2 
#
# aqua format:
# def LEU * HBX 0 1HB 2HB
#
# The constituent atoms are in PDB format (if pdb==on, the default), 
# or in standard internal format (if pdb=off).
# The pseudo atom is in the format determined by 'type'.
#

BEGIN {
    if ( type == "" ) type = "biosym"
    type = tolower( type )
    if ( type != "aqua" ) {
	print "check the output routines against mdf2AtomLIB"
	exit( 1 )
    }
    nopdb = 0
    if ( pdb != "" ) nopdb = 1
    if ( pdb == "on" ) nopdb = 0
    if ( pdb == "off" ) nopdb = 1

    "/bin/date" | getline date
    close "/bin/date"
    "printenv AQUAVERS" | getline vers
    close "printenv AQUAVERS"
    tempfil = "mdf2.tmp"
    pscount = 0
    nskip = 0
}


function setpdb( resnam, atnam,    seq, at ) {
    pdb = atnam
    if ( atnam == "HN" ) 
	return ( "H" )
    if ( match( atnam, /^H/ ) ) {
	if ( match( atnam, /[123]$/ ) ) {
	    seq = substr( atnam, RSTART, 1 )
	    at = substr( atnam, 1, RSTART-1 )
	    if ( at != "HB" ) {
		if ( match( resnam, /TYR/ ) || match( resnam, /PHE/ ) || 
		     match( resnam, /HIS/ ) || match( resnam, /TRP/ ) )
		    return ( pdb )
		if ( at == "HG" && match( resnam, /THR/ ) )
		    return ( pdb )
	    }
	    pdb = sprintf( "%s%s", seq, at )
	}
	if ( match( atnam, /X$/ ) ) {
	    at = substr( atnam, 2, RSTART-2 )
	    pdb = sprintf( "Q%s", at )
	}
    }
    if ( pdb == "QG" )
	if ( match( resnam, /VAL/ ) )
	    pdb = "QQG"
    if ( pdb == "QD" )
	if ( match( resnam, /LEU/ ) )
	    pdb = "QQD"
    return ( pdb )
}

function setiupac( resnam, atnam,    seq, at ) {
    iupac = atnam
# Code only those which change from Insight's mdf to IUPAC.
#
# Assumes only standard residues are present, look into the code to check
#  wether this will hold for your fancy AA's. OK, ACE of xplor included
# Simple to do it's ONLY CONVERTING 1/2 NOMENCLATURE TO 2/3.
#
# Protons only:
    if ( match( atnam, /^H/) && ( ! (match( atnam, /X$/ )) ) ) {
#
# alfa's (only gly not ACE of xplor)
     if (! ( match( resnam, /ACE/ ) ) ) {
      if ( match( atnam, /HA1/ ) ) return ( "HA2" )
      if ( match( atnam, /HA2/ ) ) return ( "HA3" )
     }
# beta's
     if ( match( atnam, /HB/ ) ) {
      if (!( ( match( resnam, /GLY/ ) || match( resnam, /ALA/ ) || 
               match( resnam, /ILE/ ) || match( resnam, /THR/ )  ) ) ) {
	if ( match( atnam, /HB1/ ) ) return ( "HB2" )
	if ( match( atnam, /HB2/ ) ) return ( "HB3" )
      }
     } 
# gamma's
     if ( match( atnam, /HG/ ) ) {
      if ( match( resnam, /PRO/ ) || match( resnam, /ARG/ ) || 
            match( resnam, /ASN/ ) || match( resnam, /ASP/ ) || 
            match( resnam, /CYS/ ) || match( resnam, /GLN/ ) || 
            match( resnam, /GLU/ ) || match( resnam, /ILE/ ) || 
            match( resnam, /LYS/ ) || match( resnam, /MET/ ) ) {
        
	if ( atnam == "HG1" ) return ( "HG2" ) 
	if ( atnam == "HG2" ) return ( "HG3" )
	if ( atnam == "HG11" ) return ( "HG12" ) # For Ile case
	if ( atnam == "HG12" ) return ( "HG13" )
      }
     }
# delta's
     if ( match( atnam, /HD/ ) ) {
      if ( match( resnam, /PRO/ ) || match( resnam, /ARG/ ) || 
            match( resnam, /GLN/ ) || match( resnam, /GLU/ ) ||
            match( resnam, /LYS/ ) )                               {
	if ( match( atnam, /HD1/ ) ) return ( "HD2" )
	if ( match( atnam, /HD2/ ) ) return ( "HD3" )
      }
     }
# epsilon's
     if ( match( atnam, /HE/ ) ) {
      if ( match( resnam, /LYS/ ) )                                 {
	if ( match( atnam, /HE1/ ) ) return ( "HE2" )
	if ( match( atnam, /HE2/ ) ) return ( "HE3" )
      }
     }
# if not then just return the proton.
     return ( iupac )
    }

#
# Pseudo's only. Always H?X -> Q? in case of methylenes
#                    so just do the methyls and ILE QG and a few others
    if ( match( atnam, /^H/ ) && match ( atnam, /X$/ ) ) {
# alfa's (no exceptions)
# beta's 
     if ( match( atnam, /HB/ ) && match( resnam, /ALA/ ) ) return ( "MB" )
# gamma's
     if ( match( atnam, /HG/ ) ) {
      if ( match( resnam, /VAL/ ) ) {
	if ( match( atnam, /HG1/ ) ) return ( "MG1" )
	if ( match( atnam, /HG2/ ) ) { return ( "MG2" ) 
        } else                         return ( "QG" )
      }
      if ( match( resnam, /ILE/ ) ) {
	if ( match( atnam, /HG1/ ) ) return ( "QG" )
	if ( match( atnam, /HG2/ ) ) return ( "MG" )
      }
      if ( match( resnam, /THR/ ) ) {
	if ( match( atnam, /HG2/ ) ) return ( "MG" )
      }
     }
# delta's
     if ( match( atnam, /HD/ ) ) {
      if ( match( resnam, /LEU/ ) ) {
	if ( match( atnam, /HD1/ ) ) return ( "MD1" )
	if ( match( atnam, /HD2/ ) ) { 
            return ( "MD2" ) 
	} else                         return ( "QD" )
      }
      if ( match( resnam, /ILE/ ) )    return ( "MD" )
      if ( match( resnam, /ASN/ ) )    return ( "QD" )
     }
# epsilon's
     if ( match( atnam, /HE/ ) ) {
      if ( match( resnam, /MET/ ) ) return ( "ME" )
      if ( match( resnam, /GLN/ ) ) return ( "QE" )
     }
# zeta's
     if ( match( atnam, /HHX/ ) && match( resnam, /ARG/ ) ) return ( "QZ" )
# what remains: Always H?X -> Q? this covers methylenes, arg amino, rings
     dummy = match( atnam, /X$/ ) # Get the RSTART value
     at = substr( atnam, 2, RSTART-2 )
     iupac = sprintf( "Q%s", at )
     return ( iupac )
    }
# heavy atoms are all the same so just return it.
    return ( iupac )
}

function setdiana( resnam, atnam,    seq, at ) {
    diana = setdisgeo( resnam, atnam )
    return ( diana )
}
function setdisgeo( resnam, atnam,    seq, at ) {
    disgeo = atnam
    if ( atnam == "HN" ) 
	return ( "H" )
    if ( match( atnam, /^H/ ) ) {
	if ( match( atnam, /X$/ ) ) {
	    at = substr( atnam, 2, RSTART-2 )
	    disgeo = sprintf( "Q%s", at )
	}
    }
    if ( disgeo == "QG" )
	if ( match( resnam, /VAL/ ) )
	    disgeo = "QQG"
    if ( disgeo == "QD" )
	if ( match( resnam, /LEU/ ) )
	    disgeo = "QQD"
    return ( disgeo )
}

function setxplor( resnam, atnam,    seq, at ) {
    xplor = atnam
    if ( atnam == "HRX" ) return ( "skip" )
    if ( match( xplor, /^H/ ) )
	sub( /X$/, "#", xplor ) 
#    if ( xplor == "HG#" )
#	if ( match( resnam, /VAL/ ) )
#	    xplor = "HG*"
#    if ( xplor == "HD#" )
#	if ( match( resnam, /LEU/ ) )
#	    xplor = "HD*"
    return ( xplor )
}

/@molecule/ {
    mol = $2":"
}

/@list pseudoatom/ {
    sub( mol, "" )
# remove the residue number
    sub( /_[0123456789]*:/, " " )
# check whether this pseudo atom occurred before
# and skip some special cases
    for ( i=1; i<=pscount; i++ ) {
	if ( ps[i] == $3 $4)
	    next
    }
    ps[++pscount] = $3 $4
# make an r.e. for the resid, i.e. residue name plus residue number
    resid = $3"_[0123456789]*:"
    resnam = $3
# form the name of the pseudo atom
    pseudo = $4
    if ( type == "disgeo" ) pseudo = setdisgeo( resnam, pseudo )
    if ( type == "diana" ) pseudo = setdisgeo( resnam, pseudo )
    if ( type == "xplor" )  pseudo = setxplor( resnam, pseudo )
    if ( type == "aqua" )   pseudo = setiupac( resnam, pseudo )
# write the pseudo atom in Aqua format
    if ( pseudo != "skip" )
	printf ( "def %-5s * %-5s 0 ", resnam, pseudo ) > tempfil
    else
	nskip++
# there might be a + in the residue name, so modify to make a valid pattern
    sub( /\+/, "\\+", resid )
# process the next lines, upto a blank or empty line
    while ( 1 ) {
	sts = getline
	if ( sts < 1 ) 
# error or end of file
	    exit
	if ( NF == 0 ) {
# end of definition
	    if ( pseudo != "skip" )
		printf ( "\n" ) > tempfil
	    next
	}
	if ( pseudo == "skip" )
	    continue
	gsub( resid, "" )
# write the constituent atom in standard or pdb format
	for ( i=1; i<=NF; i++ ) {
	    if ( nopdb ) 
		printf ( "%-5s ", setiupac( resnam, $i ) ) > tempfil
	    else 
		printf ( "%-5s ", setpdb( resnam, $i ) ) > tempfil
	}
    }
}

END {
# print the header
    print "$ Aqua Version", vers, "DATA: Pseudo atom definitions " toupper(type) " format"
    print "$ FILE: PseudoLIB CREATED:", date
    print "$ PSEUDO_COUNT:", pscount - nskip
    print "$ BEGIN_DATA:"
# sort the entries, and add them
    system( "/bin/sort "tempfil "; /bin/rm "tempfil )
}
